Quick tutorial - ggplot2 themes

Ever notice the gray background on ggplot2? Well that’s referred to as the default theme. There are other themes that can be called every time you createa a ggplot. Let’s take a look:

The black & white theme

library(ggplot2)

acadia = read.csv("/Users/james/Documents/Github/geog473-673/datasets/acadia.csv")

ggplot(acadia, aes(x=year, y=visitors)) +
    geom_line(col="yellow", size=5) + 
    geom_point(col="steelblue", size=3) + 
    geom_smooth(method="lm", col="firebrick") +
    labs(title="Acadia National Park Attendance", subtitle="Total Visitors per year", y="Visitors", x="Year", caption="National Park Database") +
    theme_bw()

The dark theme

library(ggplot2)

acadia = read.csv("/Users/james/Documents/Github/geog473-673/datasets/acadia.csv")

ggplot(acadia, aes(x=year, y=visitors)) +
    geom_line(col="yellow", size=5) + 
    geom_point(col="steelblue", size=3) + 
    geom_smooth(method="lm", col="firebrick") +
    labs(title="Acadia National Park Attendance", subtitle="Total Visitors per year", y="Visitors", x="Year", caption="National Park Database") +
    theme_dark()

The light theme

library(ggplot2)

acadia = read.csv("/Users/james/Documents/Github/geog473-673/datasets/acadia.csv")

ggplot(acadia, aes(x=year, y=visitors)) +
    geom_line(col="yellow", size=5) + 
    geom_point(col="steelblue", size=3) + 
    geom_smooth(method="lm", col="firebrick") +
    labs(title="Acadia National Park Attendance", subtitle="Total Visitors per year", y="Visitors", x="Year", caption="National Park Database") +
    theme_light()

The classic theme

library(ggplot2)

acadia = read.csv("/Users/james/Documents/Github/geog473-673/datasets/acadia.csv")

ggplot(acadia, aes(x=year, y=visitors)) +
    geom_line(col="yellow", size=5) + 
    geom_point(col="steelblue", size=3) + 
    geom_smooth(method="lm", col="firebrick") +
    labs(title="Acadia National Park Attendance", subtitle="Total Visitors per year", y="Visitors",x="Year", caption="National Park Database") +
    theme_classic()

The void theme

library(ggplot2)

acadia = read.csv("/Users/james/Documents/Github/geog473-673/datasets/acadia.csv")

ggplot(acadia, aes(x=year, y=visitors)) +
    geom_line(col="yellow", size=5) + 
    geom_point(col="steelblue", size=3) + 
    geom_smooth(method="lm", col="firebrick") +
    labs(title="Acadia National Park Attendance", subtitle="Total Visitors per year", y="Visitors",x="Year", caption="National Park Database") +
    theme_void()

Advanced Plotting - Synthesizing Shapefiles

Often times, the spatial data that we want aren’t in a nice, gridded lat/lon. The oblate spheroid that we call home presents some challenges when it comes to displaying spatial datasets. Often times, datasets like these are stored as separate entities - shapefiles & data tables. Fortunately, we have multiple avenues for working with this data. Let’s check out a coastal flooding product that projects potential water levels above normal for each day. [coast-flood.udel.edu]

#Load packages
library(RColorBrewer)
library(rgdal)

## Loading required package: sp

## rgdal: version: 1.4-6, (SVN revision 841)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 2.4.2, released 2019/06/28
##  Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/3.6/Resources/library/rgdal/gdal
##  GDAL binary built with GEOS: FALSE 
##  Loaded PROJ.4 runtime: Rel. 5.2.0, September 15th, 2018, [PJ_VERSION: 520]
##  Path to PROJ.4 shared files: /Library/Frameworks/R.framework/Versions/3.6/Resources/library/rgdal/proj
##  Linking to sp version: 1.3-1

library(sp)
library(ggplot2)
library(ggmap)

## Google's Terms of Service: https://cloud.google.com/maps-platform/terms/.

## Please cite ggmap if you use it! See citation("ggmap") for details.

library(scales)
library(viridis)

## Loading required package: viridisLite

## 
## Attaching package: 'viridis'

## The following object is masked from 'package:scales':
## 
##     viridis_pal

# Use the readOGR function to open a shapefile & constituents
coast.shp <- readOGR("/Users/james/Documents/Github/geog473-673/datasets/cfms_shapefiles/cfms_watersheds.shp")

## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/james/Documents/Github/geog473-673/datasets/cfms_shapefiles/cfms_watersheds.shp", layer: "cfms_watersheds"
## with 24 features
## It has 15 fields

class(coast.shp)

## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"

coast.shp@proj4string

## CRS arguments:
##  +proj=tmerc +lat_0=38 +lon_0=-75.41666666666667 +k=0.999995
## +x_0=200000 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m
## +no_defs

This projection is unique. It’s a transverse mercator with some specific lat_0 and lon_0 starting points. Notice the class of the shapefile - the underlying package controlling it’s translation to R is the sp package - A package providing classes and methods for spatial data: points, lines, polygons and grids.

NOTE: Everytime you read in a shapefile like cfms_watersheds.shp, you MUST have the .shx, .prj, .dbf, etc. all within the SAME folder. Even though we are only keying an “opening” of the .shp file, readOGR is opening them all.

ANOTHER NOTE: Notice the *@* symbol in coast.shp@proj4string. This is you query metadata associated with shapefiles.

Back to plotting…

plot(coast.shp)

# Open the dataset that corresponds to the water levels within the shapefile boxes 
coast.data <- read.csv("/Users/james/Documents/Github/geog473-673/datasets/cfms_shapefiles/water_levels.csv")
head(coast.data)

##    station       lon      lat mhhwtomsl navd88tomsl mllwtomsl     mhhw
## 1 DBOFS001 -75.14990 38.79219  2.488845   0.3779528   -2.1168 2.110892
## 2 DBOFS001 -75.14990 38.79219  2.488845   0.3779528   -2.1168 2.110892
## 3 DBOFS001 -75.14990 38.79219  2.488845   0.3779528   -2.1168 2.110892
## 4 DBOFS002 -75.30239 38.94447  2.752297   0.2916667   -2.3451 2.460630
## 5 DBOFS003 -75.30879 38.95471  2.771982   0.2877297   -2.3543 2.484252
## 6 DBOFS003 -75.30879 38.95471  2.771982   0.2877297   -2.3543 2.484252
##      mllw  maxpred maxovermhhw
## 1 -2.4948 3.622047    1.510047
## 2 -2.4948 3.622047    1.510047
## 3 -2.4948 3.622047    1.510047
## 4 -2.6368 4.048557    1.587557
## 5 -2.6421 4.048557    1.564557
## 6 -2.6421 4.048557    1.564557

# now let's find the matching key - in this case, the matching key is the "station".
coast.shp$station

##  [1] DBOFS015 DBOFS016 DBOFS003 DBOFS012 DBOFS013 DBOFS004 DBOFS007
##  [8] DBOFS009 DBOFS001 DBOFS017 DBOFS018 DBOFS022 DBOFS006 DBOFS005
## [15] DBOFS010 DBOFS011 DBOFS014 DBOFS021 DBOFS008 DBOFS001 DBOFS001
## [22] DBOFS002 DBOFS020 DBOFS003
## 21 Levels: DBOFS001 DBOFS002 DBOFS003 DBOFS004 DBOFS005 ... DBOFS022

coast.data$station

##  [1] DBOFS001 DBOFS001 DBOFS001 DBOFS002 DBOFS003 DBOFS003 DBOFS004
##  [8] DBOFS005 DBOFS006 DBOFS007 DBOFS008 DBOFS009 DBOFS010 DBOFS011
## [15] DBOFS012 DBOFS013 DBOFS014 DBOFS015 DBOFS016 DBOFS017 DBOFS018
## [22] DBOFS020 DBOFS021 DBOFS022
## 21 Levels: DBOFS001 DBOFS002 DBOFS003 DBOFS004 DBOFS005 ... DBOFS022

# notice the difference above - let's reorder the data from the shapefile and the csv data by the station, otherwise
# merging will NOT work. 
coast.shp  = coast.shp[order(as.vector(coast.shp$station)),]
coast.data  = coast.data[order(as.vector(coast.data$station)),]

# merge together the shapefile data and the csv data based on common variable, here called 'station' - we MUST use the duplicateGeoms argument because there are multiple values for the same station name
merged.coast = sp::merge(coast.shp,coast.data,by='station', duplicateGeoms = TRUE)
class(merged.coast)

## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"

# remember that package::function is another way to call a specific function. I did this becuase other packages have 'merge' functions
# and I didn't want them to be confused. 

# Now let's use brewer.pal from RColorBrewer to create a blues color pallette
mycolours <- brewer.pal(8, "Blues")
mybreaks <- seq(0,6,0.5) 
# the data we're concerned with here is shortnamed maxpred - maximum predicted water level for the next 24 hours
cut(merged.coast$maxpred, mybreaks)

##  [1] (3.5,4] (3.5,4] (3.5,4] (3.5,4] (3.5,4] (3.5,4] (3.5,4] (3.5,4]
##  [9] (3.5,4] (4,4.5] (4,4.5] (4,4.5] (4,4.5] (4,4.5] (4,4.5] (4,4.5]
## [17] (4,4.5] (4.5,5] (4.5,5] (4.5,5] (4.5,5] (4.5,5] (4.5,5] (4.5,5]
## [25] (4,4.5] (4.5,5] (4.5,5] (4,4.5] (4,4.5] (3.5,4] (4,4.5] (4,4.5]
## 12 Levels: (0,0.5] (0.5,1] (1,1.5] (1.5,2] (2,2.5] (2.5,3] ... (5.5,6]

mycolourscheme <- mycolours[findInterval(merged.coast$maxpred, vec = mybreaks)]

The data is prepared, now we just need to plug in everything to plotting functions. First up is spplot from the sp package:

spplot(merged.coast, "maxpred", par.settings = list(axis.line = list(col ="transparent")), main = "Projected Water Levels (Feet)", cuts = 5, col ="transparent", col.regions = mycolours)

Now we use the tmap function - very similar to ggplot2. tmap == Thematic Map Visualization. Thematic maps are geographical maps in which spatial data distributions are visualized. This package offers a flexible, layer-based, and easy to use approach to create thematic maps, such as choropleths and bubble maps. It is based on the grammar of graphics, and resembles the syntax of ggplot2.

library(tmap)
tm_shape(merged.coast) + 
  tm_polygons(col='maxpred', title = "Projected Water Levels", palette = "Spectral") + tm_style("classic") + tm_scale_bar(position = c("left", "bottom"))

Ok, so those are cool, but what about using ggplot2? In order to do this, we’ll need to use geom_sf from ggplot2. This enables us to plot sf objects. sf is a package for Simple Features. It’s one of the other big players for spatial data wrangling in R that a lot of other, fancier packages that we use in this class are built on. We’ll need to load in sf before we can convert our shapefile into an sf object.

library(sf)

## Linking to GEOS 3.7.2, GDAL 2.4.2, PROJ 5.2.0

states <- readOGR("/Users/james/Documents/Github/geog473-673/datasets/ne_10m_admin_1_states_provinces/ne_10m_admin_1_states_provinces.shp")

## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/james/Documents/Github/geog473-673/datasets/ne_10m_admin_1_states_provinces/ne_10m_admin_1_states_provinces.shp", layer: "ne_10m_admin_1_states_provinces"
## with 4594 features
## It has 83 fields
## Integer64 fields read as strings:  ne_id

class(states)

## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"

states = st_as_sf(states)
class(states)

## [1] "sf"         "data.frame"

# convert it to an sf object
gg_merged = st_as_sf(merged.coast)
class(gg_merged)

## [1] "sf"         "data.frame"

# now we use geom_sf since our merged.coast is now a simple feature
ggplot(data = states) + geom_sf() + theme_void() +
  geom_sf(data = gg_merged, aes(fill = maxpred)) +
  coord_sf(xlim = c(-76.1, -74.7), ylim = c(38.3,40), expand = FALSE) +
  scale_fill_distiller(palette = "Blues", direction= 1)

Looks pretty good even though the shapefile isn’t quite as sharp as we’d like. Besides the shapefile which we added to the plot via geom_sf(), take notice of the theme_void()…this plot will NOT WORK without theme_void(). There is ongoing github discussions with some of the developers at ggplot2 to figure out why, but for now just make sure you use theme_void to add shapefiles to the same plot in R and have them work properly. Also take note of the scale_fill_distiller function. This is an easy way to throw in a RColorbrewer color pallette. In this case we used the Blues colorpallete since we’re dealing with water but you can make it any RColorbrewer pallette you want. Also notice the direction=1. This reverses the order of the blues color pallette. I can’t explain why it’s not -1…Sidenote - you can also use scale_fill_viridis() which has some different default pallettes such as - “magma”, “plasma”, and “inferno”.

Mapview

Saving the easiest for last…mapview is an awesome package that I didn’t find until recently. As always though, a house with a foundation is more sound than one without. It’s good to know about the other packages above!

library(mapview)
mapview(merged.coast['maxpred'], col.regions = mycolours)

Saving a new shapefile

So we combined a csv with a shapefile to make that plot. We don’t need to duplicate that hard work again because we can actually save our merged shapefile/csv as an ENVI shapefile.

writeOGR(obj, dsn, layer, driver, dataset_options = NULL, layer_options=NULL, verbose = FALSE, check_exists=NULL, overwrite_layer=FALSE, delete_dsn=FALSE, morphToESRI=NULL, encoding=NULL, shp_edge_case_fix=FALSE)

writeOGR(obj = merged.coast, dsn = "/Users/james/Downloads/coast_files/merged.coast", layer = "coast-rgdal", driver = "ESRI Shapefile")

Week 5 Assignment - Replicate these Plots

Using the ne_10m_parks_and_protected_lands shapefile dataset, map the names of the protected lands in the Pacific Northwest using the methods listed below. You will need to download all of the contents of the ne_10m_parks_and_protected_lands shapefile folder.

spplot
tmap
ggplot2
mapview

Create the 4 plots above and make them look as close as possible to the ones below.

## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/james/Documents/Github/geog473-673/datasets/ne_10m_parks_and_protected_lands/ne_10m_parks_and_protected_lands_area.shp", layer: "ne_10m_parks_and_protected_lands_area"
## with 61 features
## It has 8 fields
## Integer64 fields read as strings:  scalerank

## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/james/Documents/Github/geog473-673/datasets/ne_10m_admin_1_states_provinces/ne_10m_admin_1_states_provinces.shp", layer: "ne_10m_admin_1_states_provinces"
## with 4594 features
## It has 83 fields
## Integer64 fields read as strings:  ne_id

## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"

## [1] "sf"         "data.frame"

Extra Credit - 3 points !!!

Using the ne_10m_parks_and_protected_lands shapefile dataset, create the following…

spplot with statelines and new color scheme
tmap with statelines and new color scheme
ggplot2 with legend containing only protected lands in the Pacific NW
mapview with new color scheme

Submit resulting plots and code as “ec_plottype.png” with the above assignment

Week5 - Shapefiles